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Abstract 

The four-body bound state with two-body interactions is formulated in Three-Dimensional ap- 
proach, a recently developed momentum space representation which greatly simplifies the nu- 
merical calculations of few-body systems without performing the partial wave decomposition. 
The obtained three-dimensional Faddeev-Yakubovsky integral equations are solved with two-body 
spin-independent and spin-averaged potentials. This is the first step toward the calculations of 
four-nucleon bound state problem in Three-Dimensional approach. Results for four-body binding 
energies are in good agreement with achievements of the other methods. 
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I. INTRODUCTION 



The bound state of few-body systems seems to be an ideal laboratory to determine two-, 
three- and four-body nuclear forces. The studies of the four-body bound state properties 
for the case of few-body interactions have received increasing attention theoretically and 
experimentally in recent years. Although the four-body bound state poses a challenging 
problem numerically, because of presence of fourth body, its investigation promises insights 
into the rich structure of nuclear interactions. To this aim one requires an accurate and 
reliable method to obtain the full solution of four-body bound state in a straightforward 
manner. 

The four-body bound state calculations are carried out by different methods to solve 
the nonrelativistic Schrodinger equation such as, the Coupled-Rearrangement-Channel 



Gaussian-basis Variational(CRCGV) 



Gaussians 



fl-fl, 



the Stochastic Variation(SV) with correlated 
12], the Hyperspherical Harmonic variational (HH) [3]- [3], the Green's Func- 



tion Monte Carlo(GFMC) (l9|-[22j, the No-Core Shell Model(NCSM) |23|-|28j, Effective 
Interaction Hyperspherical Harmonic(EIHH) 29],j3o| and the Faddeev-Yakubovsky (F-Y). 
In the last method the nonrelativistic schrodinger equation is transformed to two coupled 
sets of finite number of coupled equations in three variables for the F-Y amplitudes. The 



3l| - 34( and in momentum 



calculations based on F-Y are performed in configuration space 
space 35J- 43j after a partial wave (PW) expansion, where the algebraic and algorithmic 



steps can be quiet involved. Though a few partial waves often provide qualitative insight, 
modern four-body calculations need 1572 or more different spin, isospin and angular mo- 



mentum combinations 



42], 43j. It appears therefore natural to avoid a PW representation 



completely and work directly with vector variables. On this basis in recent years W. Glockle 
and collaborators have introduced the Three-Dimensional (3D) approach which greatly sim- 
plifies the two- and three-body scattering and bound state calculations without using PW 
decomposition [44j-[53[]. 

Our aim in this paper is to extend this approach for four-body bound state with two- 
body interactions, we work directly with vector variables in the F-Y scheme in momentum 
space. Here the dependence on momentum vectors, i.e. the magnitudes of momenta and the 
angles between the momentum vectors, shows that the full solution can be reached exactly 
and simply whereas the PW representation of the amplitudes leads to rather complicated 
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expressions 



541 ] . The calculations of four-body bound state with the three-body interactions 



is currently underway and it will be reported elsewhere. As a simplification we neglect spin 
and isospin degrees of freedom here and study the four-boson bound state. So this work 
is the first step in the direction of solving the four-nucleon bound state problem without 
performing the PW decomposition. 

Recently the four-boson bound state has been studied with short-range forces and large 
scattering length at leading order in an Effective Field Theory approach j^j-js^], but this 
investigation is also based on PW decomposition and the interactions are restricted to only 
S-wave sector. 

This paper is organized as follows. Section [TT] reviews the F-Y equations for four-boson 
bound state. In section II I II we represent the coupled F-Y equations as function of momen- 
tum vectors. In section IIVI we discuss our choice for independent variables of momentum 
and angle variables for the unknown amplitudes in the equations and in their kernels, where 
this new representation(3D) is contrasted with traditional PW representation. Section [V] de- 
scribes details of our algorithm for solving coupled F-Y three-dimensional integral equations. 
In section |VT] we compare our results for three- and four-boson binding energies to results 
obtained from other techniques. In order to test our calculation we investigate the stability 
of the eigenvalue of the Yakubovsky kernel with respect to the number of grid points and we 
calculate the expectation value of Hamiltonian operator. Finally we summarize in section 
IVHI and provide an outlook. 



11. FOUR-BODY BOUND STATE EQUATIONS 

The bound state of four identical particles which interact via pairwise forces Vij(ij 

1 2 . 13, 14, 23, 24 and 34) is given by Schrodinger equation which reads in integral form: 

|*) = Gfe$><,-|¥) (1) 

i<j 

Here the free four-body propagator is given by Gq = (E — Hq)' 1 , and H stands for the free 
hamiltonian. Introducing Yakubovsky components \^/) = with = G^Vijl^) 

leads to the six coupled integral equations: 

\M = Got tJ l^«> ( 2 ) 
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The operator Uj describes the two-body t— matrix in the two-body subsystem ij. We can 
rewrite Eq. (jSJ) as: 

= Gotijdipik) + \ip u ) + \tfj jk ) + + \ipki)) (3) 

Among various possibilities to decompose \4>ij) into three F-Y components we choose the 
following one: 

\Ajk,l-,ij) = GoUjdipik) + \i>jk}) 

\Ajl,k;ij) = GotijQlpu) + 

\i>ij,H;ij) = G Q tij\ip k i) (4) 

The F-Y component \ifjijk,i-,ij) i\^ij,ki;ij)) belongs to a 3 + 1 (2 + 2) partition. They fulfill the 
following relation: 

= \i>ijk,l;ij) + \4>ijl,k;ij) + \AjM\ij) ( 5 ) 

The multiple indices for the F-Y components denote the two-body followed by the 3 + 1 or 

2 + 2 fragmentation. It is easily seen that every component contains two 3 + 1 type 
chains and one 2 + 2 type chain, therefore total wave function \^f) contains twelve different 

3 + 1 type chains and six 2 + 2 type chains. So altogether one has eighteen F-Y components. 
If we consider identical particles (here bosons, since we are omitting spin), the four-body 
wave function has to be totally symmetric. As a consequence all twelve components of 
3 + 1 type are identical in their functional form and only the particles are permuted. The 
same is true for the six components of 2 + 2 type. Thus it is sufficient to consider only two 
independent F-Y components corresponding to the 3 + 1 and 2 + 2 partitions, 

= I^123,4;12> 

1^2) = |^12,34:12> (6) 

After the straightforward derivation the 18 coupled F-Y components shrink to two coupled 
F-Y equations: 

|^l>=Gotl 2 P[(l + P34M> + l^>] 

M = G t Vi P[{l + P si )M + \ih)] (7) 
Where Pij is the permutation operator between the i — th and j — th particle, and 

P = P12P23 + P13P23 

P = P13P24 (8) 
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The total four-body wave function is then given as: 

|*) = (1 + P + P U P + P)[(l + P 34 )|^) + |V> 2 >] (9) 

The symmetry property of l^i) under exchange of particles 1 and 2, and \ip 2 ) under separate 
exchanges of particles 1,2 and 3,4 guarantee that \^) is totally symmetric. 

We would like to add the remark that another derivation of F-Y components is also 
possible |35j . In this representation two transition operators which follow the subcluster 
Faddeev-like equations have been introduced as a function of two-body transition operator 
ti2- Consequently the kernel of coupled Yakubovsky integral equations contains two sub- 
cluster kernels that should be evaluated by Pade technique. So its numerical calculations is 
more complicated and time consuming in comparison to above derivation. 



III. MOMENTUM SPACE REPRESENTATION OF FADDEEV- YAKUBOVSKY 
EQUATIONS 

In order to solve the coupled equations (7), in momentum space we introduce standard 
Jacobi momenta sets corresponding to both 3 + 1 (123, 4; 12) and 2 + 2 (12, 34; 12) chains: 



h - k 2 



2 r fci + fo 

u 2 = -(fc--^— ) ^ 

3 r hi + k 2 + k 3 > 
u 3 = ^ , 

h - k 2 



V2 



2 

ki + k 2 k 3 + k 4 
2 2 



k 3 — k A 

v 3 = 2 ^ 10 ) 

Then we introduce the four-body basis states corresponding to each Jacobi momenta set: 

\ut u 2 u 3 ) 

W1V2V3) (11) 
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Both basis states are complete in the four-body Hilbert space: 

J D 3 A \A 1 A 2 A 3 ){A X A 2 A 3 \ = 1 (12) 

Where A4 indicates each one of Ui and Vi vectors and D 3 A = d 3 A\ d 3 A2 d 3 A 3 . Also they are 
normalized according to: 

(A A 2 A 3 \A[ A' 2 4) = 5 3 {A 1 - A[) 5 3 (A 2 - A! 2 ) 5 3 (A 3 - A' 3 ) (13) 

Clearly the basis states \uiu 2 u 3 ) are adequate to expand F-Y component [tpi) and corre- 
spondingly the basis states \vi v 2 v 3 ) are adequate for |^ 2 ). Let us now represent coupled 
equations, Eq. ([7]), with respect to the basis states have been introduced in Eq. ( fTTT) : 



(^U.UslA) = J W , (M 1 ii 2 ii3|G'otP(l + P34)|^l«2W30(«l«2«3lV'l) 

+ / D 3 v' (u 1 u 2 u 3 \G tP\v[v' 2 v' 3 )(v[v' 2 v' 3 \tp 2 ) 



(viv 2 v 3 \i{j 2 } = j D 3 u' (viV 2 v 3 \G Q tP{l + P 34t )\u[u 2 u 3 )(u[u 2 u' 3 \ip 1 ) 

+ /DV^^^IGotPl^^^)^^^!^) (14) 



It is convenient to insert again the completeness relations between permutation operators, 
it results: 



(uiu 2 u 3 \ijji) = I D 3 u ' j D 3 u'' (uiu 2 u 3 \G tP\u[u' 2 u' 3 } 

x (u[ u' 2 u' 3 \{l + P u )\u'l ul u'i) (u'i u'i u'i^x) 
D 3 u J £>V (utfauslGotPlu'iU^u'z) 
x (u'i u 2 u 3 \v[ v 2 v' 3 ) (v[ v 2 v' 3 \if; 2 } 



+ 



(v!V 2 v 3 \^ 2 ) = j D 3 v' j D^u'iv^VslGotPlv^v'^) 

x(v[ v 2 v 3 \(l + P u )\u[u 2 u 3 )(u' 1 u 2 u 3 \4> 1 ) 
+ /DV^^^IGotPl^^^^W^al^) (15) 
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For evaluating the coupled equations, Eq. (|T5|) . we need to evaluate the following matrix 
elements: 



{Hi u 2 u^GotPlu^ u 2 u 3 ) (16) 

(v 1 v 2 v 3 \G tP\v' 1 v' 2 v 3 ) (17) 

(u\u' 2 u' 3 \(l + P 3A )\u^ul) (18) 

(v[v' 2 v' 3 \(l + P 3i )\u[u' 2 u' 3 ) (19) 



For evaluating the first term, Eq. (1161) . we should insert again a completeness relation 
between the two-body t— matrix operator and permutation operator P as: 

(mi u 2 u^GotP^ u 2 u 3 ) = ul —z ^ 

m 4m 3m 

x J P) 3 u" (u\ u 2 u 3 \t\u'^ u 2 u 3 ) (u'l u 2 u' 3 \P\u'i u 2 u' 3 ) (20) 

Where the matrix elements of two-body t— matrix and permutation operator P are evaluated 
separately as: 

(iZi u 2 u 3 \t\u'[ u 2 u 3 ) = 6 3 (u 2 - u 2 ) 5 3 (u 3 - u" 3 ) (ui|t(e)|u'/) 

e = E — — ^ ^ 21 

4m 3m 

(MiWaWglPlMiUaUg) = <5 (u 3 -u 3 ) 

x { 5 3 « + ^ - ^) 8\u' 2 > + u[ + \u' 2 ) 

+ Pffl + \u' x + \u' 2 ) 5 3 (u' 2 > - u[ + \u' 2 ) } 

(22) 

For evaluation the matrix elements of permutation operator P we have used the rela- 
tion between Jacobi momenta in different two-body subsystems (312,4; 12), (231,4; 12) and 
(123,4; 12). Inserting Eqs. fl2B) and ([22]) into Eq. ([20]) leads to: 

5 3 {u 3 - u' 3 ) 



(ux u 2 u 3 \G tP\u' 1 u' 2 u' 3 ) 



m 4m 3m 

{ 5\u 2 + u[ + -ti' 2 ) {u x \t{e)\-u 2 + u' 2 ) 
+5 3 {u 2 - u[ + -u' 2 ) {ux\t{e)\^-u 2 - u' 2 )} 



(23) 
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Representation of the second term, Eq. (fl7|) . follows the similar steps: 

(M^\G tPWW 3 ) = — ^ / D 3 v"(v 1 v 2 v 3 \t\v'(v' 2 'v^ 

JTj "l v 2 "3 J 

m 2m m 

x fflWn\PK*W (24) 
The matrix elements of two-body t— matrix and permutation operator P are evaluated as: 

{vxV 2 vMv'[v'ivl) = 5\v 2 -vl)6\v,-v^){v 1 \t{e*)\v'[) 

e * = E - ii - V l (25) 
2m m 



(t'lvU^PWW's) = S 3 (v';-v' 3 )6 3 (^ + v' 2 )6 3 (^-v' 1 ) (26) 
Inserting Eqs. ([25]) and ([26]) into Eq. ([24]) leads to: 

(VMlGatPMVM) = * 3 fo + *2)* 3 fo-g'i) (^1^)1^) (27) 

E v i "2 ^3 

m 2m m 

For evaluation the third term, Eq. ( fl8l) . we should use the relation between Jacobi momenta 
in different chains (123,4; 12) and (124,3; 12), which leads to: 

(u[u' 2 u' 3 \(l + P M )\u';u'>u'^ = 5 3 {u' x -u'{) 

x {8 3 (u 2 -u'£6 3 (u 3 -u'£ 

+5 3 (u> 2 - \u'i - % -u'i) 5 3 (u> 3 - ul + \ul) } 



(28) 

And finally for evaluation the fourth term, Eq. ( fT9j) . we should use the relation between 
Jacobi momenta in two naturally different chains (123, 4; 12) and (12, 34; 12), which leads 
to: 

(*W' a \(i + p M )\iil*' a ii' a ) = aVi-tfi) 

x { 5 3 {u' 2 + -v 2 - ~v' 3 ) 5 3 (u' 3 + -v' 2 +v' 3 ) 
+8 3 (u' 2 + 2 -v' 2 + 2 -v' 3 ) 6 3 {u' 3 + -v' 2 - v' 3 ) } 

(29) 
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Finally inserting Eqs. (551). (f2Tj) . (j25jl and p^j) in Eq. (fT5|) yields: 

If 1 

(wiM 2 M 3 |^l) = ^ ^2 / ^V 2 (t*l|* a (e)b«2 + «' 2 ) 

E — -J- 2. — — 3 J Z 

m 4m 3m 

X { («2 + 7^2 #2 "sIV'l) 

+ (« 2 + -« 2 + ^U 3 u' 2 - ^MslV'l) 

1 2 12 

+ (w 2 + -m 2 - u' 2 - -m 3 -w' 2 - -U 3 \^ 2 ) } 

(^^^3^2) = j^-*-^ ^ / <^ 3 (^( e *)K 7 3> 

£j _ Hi. _12_ 12. J 

m 2m m 

2 2 1 

X { 2 (#3 7^2 + -# 3 ~ 

+ (v 3 -v 2 v' 3 \ip 2 ) } (30) 

Here {a\t s (e)\b) generally represents the symmetrized two-body t— matrix which is defined 

as, 

(S\t s (e)\b) = (a\t(e)\b) + (a\t(e)\ - b) (31) 

We would like to mention that the so obtained F-Y amplitudes fulfill the below symmetry 
relations, as can be seen from Eq. fl30|) : 



(VlV 2 V 3 \ll) 2 ) = {~ViV 2 V 3 \tp2) 

(viv 2 v 3 \^ 2 } = (viv 2 -v 3 \ip 2 ) (32) 

From the F-Y components and \i/j 2 ) the four-body wave function is obtained by adding 
the components defined in different 3 + 1 and 2 + 2 type chains as given in Eq. ([9l). After 
evaluating the permutation operators P, P and P34 the wave function is given as: 

I*) = + |*2) (33) 



9 



Where (^2)) is corresponding to all 3 + 1 (2 + 2) type chains: 

(ui u 2 usl^i) = 

{ (ui u 2 Msl^i) 

— 1 _ 3 1 _ 

+ - 7^2 "l - " 3 '^ 

+(^1 + ^2 -^1-^2 « 3 hM> 123+4 

+ { {Ui ~U 2 + -« 3 M 2 - gtWl) 

— 1 _ 1 _ 2 14 1 

+ (^T^i - 7^2 - Mi - -m 2 - 77 m 3 m 2 - -M^l) 
2 4 3 o 9 3 

— 1 _ 1 _ 2 14 1 

+(^T^i + 7^2 + 7M3 - Mi - 77M2 - m 2 - -m 3 |^i) } 124+3 
2 4 3 o9 3 

+1 (g"! ~ 4^2 - -ui - -u 2 + -u 3 - mi - -m 2 - -nam) 
1 2 2 2 4 1 1 

+ ^ 2 _ 3^ 3 3" 1 ~ 3^ 2 ~ 9^ 3 ~~ ~ 2^ 2 ~ 3^ 3 ^ 

U 2^ U 5^ 4^ U U , , m 

+ + 4 W2 + -gMl+gM2--M3 ~ U! - -U 2 - -U^!) } 

+1 ( _ 2 Ml ~ I" 2 3 Ml ~~ 6 M2 9 M3 Ul ~ 2 U2 ~ s U3 ^ 
1 2 2 2 4 1 1 

+ (7^2 - -u 3 - -ui - -m 2 - -m 3 Mi - -m 2 - -m 3 |V>i) 

/I- 2 - 5 - 1- , , m 234+1 

+\2 Ul + 4 M2 + 3 W3 3 Ml + q U2 ~ g M3 Wl ~ 2 U2 ~ s U3 ^ 



12+34 



(Mi v 2 v 3 \^ 2 ) = 
{ (vi v 2 v 3 \ip 2 ) 

+ (v 3 - V 2 V!\l(j 2 ) } 

+{ (t^ 1 + 2 V2 _ 2^ 3 Vl + V3 ~~ 2 Ul + 2^ 2 + 2 V3 '^ 2 ' 
— 1 1 1 1 1 1 

+ (^-^1 + 2^ 2 + 2^ 3 - ^ 1- ^ 3 Tj^i + 2^2 - 2^l^ 2 ) } 13+24 

111 111 

+ (--^l + -^2 - -V 3 -V X + V 3 -Mi + -V 2 + -V 3 |^ 2 ) } 14+23 
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Each curly bracket contains all possible chains in the subsystem which is indicated with 
corresponding superscript. Already here we see that: 



(wiM2M 3 |^l) = (-Ui W2M 3 |*i) 

(viv 2 v 3 \^> 2 ) = (viv 2 -t? 3 |*2> (36) 



Eq. (1361) is satisfied if the F-Y components fulfill the expected symmetries in Eq. (1321) . 

IV. CHOOSING THE COORDINATE SYSTEMS 

The F-Y components \tpi(Ai A 2 A 3 )) are given as a function of Jacobi momenta vectors 
and as a solution of coupled three-dimensional integral equations, Eq. ( l30l) . Since we ignore 
spin and isospin dependencies, the both F-Y components \ipi(Ai A 2 A 3 )) are scalars and 
thus only depend on the magnitudes of Jacobi momenta and the angles between them. The 
first important step for an explicit calculation is the selection of independent variables. One 
needs six variables to uniquely specify the geometry of the three vectors A±, A 2 and A3, 
which are shown in Fig. [TJ Having in mind that with three vectors one can span 2 planes, 



i.e. the A3 — A\ plane and A3 — A 2 plane, a natural choice of independent variables is 53]: 



Ai = |A x | A 2 = \A 2 \ A 3 = |A 3 | 

Xl = A 3 X x 2 = A 3 .A 2 x\ 2 = (A 3 x A 1 ).{A 3 x A 2 ) (37) 

The last variable, x\ 2 , is the angle between the two normal vectors of the A3 — A\ plane 
and the A 3 — A 2 plane, which is explicitly related to the angle between Ai and A 2 vectors 
as: 

3 I1.I2 -(AVlxXlg.A) 

X 12 — 1 1 

sJl-^.A.ysJl-^.Atf 
Therefore in order to solve Eq. ( |30l) directly without employing PW projection, we have to 
define suitable coordinate systems. As shown in Fig. [2} for both F-Y components we choose 
the third vector parallel to Z— axis, the second vector in the X — Z plane and express the 
remaining vectors, the first as well as the integration vectors, with respect to them. We have 
the magnitudes of vectors as well as the following angle relations as variables: 
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FIG. 1: The geometry of three vectors A±, A2 and A3 relevant in the four-body bound state 
problem. The independent angle variables x\, X2 and x\ 2 as defined in Eq. (|37p are indicated. The 
dashed arrows represent the normal vectors (A3 X A\) and (A3 x A2). 





7^Y 



FIG. 2: These figures show the geometry of both vector sets (ui, U2, U3) and (v%, V2, #3) relevant 
in the F-Y componets. The third vectors U3, V3 have been chosen parallel to Z— axis, the second 
vectors U2, v% i n the X — Z plane and the first vectors u\, v\ are free in the space. The independent 
angle variables x\, X2,x\ 2 and X±, X2,Xf 2 as defined in Eqs. (|39|) and (f40l) are indicated. 
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X\ = U 3 .U\ = cos(i?i) 

X 2 = U 3 .U 2 = COs($ 2 ) 

x\ 2 = u[ y .U^ = cos(<^i) 

2/12 = U\.u 2 = X\X 2 + \Jl - x\\Jl - x\ cos(v?i) 

x' 2 = u 3 .u' 2 = cos(i9 2 ) 

t/ 12 ' = Mi.M 2 = X\X 2 + \J 1 — x\\J 1 — x' 2 2 COs((/?i — </? 2 ) 

y 22 > = m 2 .m 2 = x 2 x 2 + sjl-xl^jl- x 2 2 cos((^ 2 ) (39) 



*1 


= V 3 .V X = 


cos(6'i) 


^2 


= v 3 .v 2 = 


cos(6> 2 ) 


A 12 


^<cy ^ay 

= v 1 .v 2 


= COs(0i) 


^12 


= Vl-V 2 = 


X ± X 2 + \J 


*3 


= Vs.tfz = 


cos(^) 


^13' 


= Vi.ifs = 


x 1 x' 3 + y/: 




= V 2 -^3 = 


X 2 X' 3 + J 



(40) 



With this choice of variables the matrix elements of F-Y components are given as: 



(Hi u 2 u 3 \^i) = 4>i (ui u 2 u 3 Xi x 2 x 3 



12/ 
3 



(viv 2 v 3 \ip 2 ) = ^ 2 (v 1 v 2 v 3 X 1 X 2 X 12 ) (41) 

Furthermore (ui\t s (e)\^u 2 +u 2 ) and (vi\t s (e*)\v' 3 ) are also scalar functions, and then can be 
written in the following form: 

(ui\t s (e)\^u 2 + u' 2 ) = t s (ui,^u 2 + u' 2 ;e) = t s (wi, vr, x; e) 

(vi\t s (e*)\v' 3 ) = t a ^ 1 ,v' 3 ;e*) = t s (v 1 ,v' 3 ,Y 13:] e*) (42) 

Where 



vr = \^u 2 + u' 2 \ = y^ul + u' 2 2 + u 2 u' 2 y 22 > 

x = u 1 .(-u 2 + u' 2 ) = ^(-u 2 y 12 + u' 2 y 12 >) (43) 

Z 7T I 
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The more complex dependencies appear under the integrals in Eq. (1301) for magnitude and 
angle variables of F-Y components. According to Eq. fl37|) and Eq. fl38|) they are given as: 



(u 2 + ^u' 2 u' 2 u 3 \ipi) = ^i(tti u' 2 m 3 x 12 xi 3 x^) 



(m 2 



1 



-it , 



-U 2 H W 3 M 
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2 " gW 3 |<0l) = ^l(7Ti 7T 2 7T 3 X 22 X 23 



(«2 + 



«2 



2 1.2 

2^3 -M 2 - -W 3 |V^) = ^(TTl 7T4 7T 5 X 32 X 33 X 



7Tg "> 
7Ti7r4 y 



<u 3 + - v' 3 \iJi) = ip!(v 3 Si S 2 X 12 X 13 X v ^] 

(v 3 - v 2 v' 3 \ip 2 ) = i(j 2 (v 3 v 2 v 3 X 22 X 23 X^ V2 ) 



(44) 



Where the shifted arguments are: 



7Tl 
7T 3 
7T5 

Si 

s 2 



U 2 + 2 U 2l 



M 2 + 1 U 2 + u 2 u' 2 y 22 > 



-u , 



- Q us\ 



1 .„ 64 _ 16 , . 

g U 2 + g^3 + ^f U 2 U 3 X 2 



U 2 - -u 3 \ 



-u 2 - -u 3 \ 



«2 + q U 3 



-u 2 u 3 x 2 



-u , 



U 2 + ^ u 3 + ^U' 2 U 3 X 2 
\u 2 + - ^2«3Z2 



3^2 + -V3I = 3 \hl + v 'i + 2v 2 v' 3 Y 23 , 



2^2 -v 3 l 



-v 2 2 + ^ 2 - « 2 w 3 y 23 ' 



(45) 



x 



Z11 
£12 

Zl3 

J 7T1W2 



1 1 1 

(u 2 + -u' 2 ).u 2 = —(u 2 y 22 > + -w 2 ) 

Z TV 1 I 



(u 2 + -u' 2 ).u 3 = 
u 2 .u 3 = x 2 

X\i — Xi 2 Xi 3 



TV 1 I 



-12 



2 

13 



(46) 



14 



X 



7T1 7T4 
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£21 = (U 2 + ^U' 2 )-(^U' 2 + -U 3 ) 

1 A i 8 1 a 4 / m 

= {^u 2 u 2 y 22 > + -u 2 u 3 x 2 + -u 2 + -u 2 u 3 x 2 ) 

" r ' r 

£22 = (u 2 + -u' 2 ).{u' 2 - -u 3 ) 

1 111 

= (u 2 u' 2 y 22 > - -u 2 u 3 x 2 + -u 2 - -u' 2 u 3 x' 2 ) 

7Ti7T3 " "~ 6 A O 

£23 = (^ / 2 + -M 3 ).(M / 2 -^M 3 ) 
1 A.J , .., 8 



'o*4 + 7T«2«34 - 7^%) 



7r 2 7r 3 '3 z 9 z z 27 

7T3 £21 — ^22^23 



X 



23 



1 2 

£31 = (U 2 + -U' 2 ).(-U' 2 - -U 3 ) 

— 1 2 11 

■(u 2 u' 2 y 22 > + -u 2 u 3 x 2 + -w 2 2 + -M 2 -U3£2) 



7Ti7T4 3 2 3 

^32 = [U 2 + -U 2 ).{-U' 2 - -U 3 ) 

11 2 11 

{-u 2 u' 2 y 22 > - -u 2 u 3 x 2 + -w 2 2 - -u' 2 u 3 x 2 ) 



7Ti7T5 2 ~ 3 

' 2 1 ~ 2 
£33 = {-u' 2 ~ -u 3 ).(-u 2 - -u 3 ) 

-1 A , 1 / / 4 2^ 

: M 2 - o M 2«3£ 2 - 7T%) 



7r 4 7r 5 v 2 z 3 z z 9 

£31 — £32£33 



x 



33 



X n = M^v 2 + ^ 3 ) = ^-(v 2 X 2 + v' 3 X' 3 ) 

^12 = ^-(^ - #3) = ^(2^X2 - V' 3 X' 3 ) 

2 ' 2 1 1 - 1 1 

*13 = (7^2 + 3^3)-(2^2 - V' 3 ) = ^-(-u* - 2^2^23' 
^£2 -^"ll ~ ^12^13 
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X 2X = v 3 .(-v 2 ) = -X 2 

X 22 = ^3-^3 = X 3 

X 23 = (-v 2 ).tf 3 = -Y 23 > 

xv > 3 _ x 21 - x 22 x 23 

v 3v 2 n YT R yz~ ^ ' 

V 1 ~ A 22V 1 ~ A 23 

These considerations lead to the explicit representation for the F-Y components and 
l^2>: 



T/>l(Wi U 2 U 3 X X x 2 xf 2 ) 



m 4m 3m 



/*00 /"2"7T 

/ du' 2 u 2 / dx 2 / dip 1 , 
Jo J-i Jo 



t s (ui,Tt,x;e) 

X { ^i(7Ti w 2 U 3 X 12 Xi 3 x^ u ,J 
+lfj 1 (7T 1 7Y 2 7T 3 £ 22 X 23 X^J 
+^ 2 (VT1 7T 4 7T 5 £ 32 X33 X™ VA ) } 



i>2(vi v 2 v 3 X 1 X 2 X\ 2 ) 



E - £ - -^i - 

m 2m m 



7 / /2 

du 3 v 3 



+1 
1 



dX 3 



2- 



x { 2^(v 3 Si S 2 X 12 X 13 Xj s J 
+^2(^3 v 2 v 3 x 22 x 23 xgj } 



(51) 



The above coupled equations, Eq. (lo"Tj) . is the starting point for numerical calculations, and 
the details will be described in the next section. In a standard PW representation Eq. (ITS]) 
is replaced by two coupled sets of finite number of coupled integral equations |4Gj: 



16 



3„.W 



im ii2 n:i^i\i'0 -- 1 ,): ''"'Y1 1 ,): '""Y1 1 /; ''' 

ft/ II J III J 

Qj ot{ a,\ 

\u x U2 U3 a x \Got\u x u 2 u 3 a x ) [u x u 2 u 3 a 1 \P\u x u 2 u 3 a x ) 

I 11 11 11 li\/-t , r> m ill 111 m 11 1 \ 1 111 111 111 ill] i \ 

{u x u 2 u 3 a x \{l + F34JF1 u 2 u 3 a x ){u x u 2 u 3 a x 



E 



(vi v 2 v 3 a 2 \ip 2 ) 



(u x u 2 u 3 a x \Got\vl x u 2 u' 3 a' x ) (u' x u 2 u 3 a' x \P\u'[ u 2 u 3 a") 

I 11 11 11 ill I 1 1 1 \ I 1 I I I I ; \ 

{u x u 2 u 3 a x \v x v 2 v 3 a 2 ) {v x v 2 v 3 a 2 \ip 2 ) 

E / flV E / BV 'E / dV 

/ rj I/O It/ 

a 2 a 2 a l 

I j.1 / / / I \ I I I I 1 \ r>\ 11 11 11 n\ 

{v x v 2 v 3 a 2 |G *K v 2 v 3 a 2 ) {v x v 2 v 3 a 2 \P\v x v 2 v 3 a 2 ) 

I 11 11 11 111 1 1 1 1 \ 1 1 1 1 1 1 / 1 i r> \ i \ 

{v x v 2 v 3 a 2 \u x u 2 u 3 a x ){u x u 2 u 3 a x \(l + P^m) 
+ J2J DV 5Z j D 3 v" (v x v 2 v 3 a 2 \G t\v' x v' 2 v' 3 a' 2 ) 

a 2 a 2 

(v' x v' 2 v' 3 a' 2 \P\v'{ v'i v'i a® « v' 2 ' v' 3 > oflfo) (52) 

Where a x = (hl2)h, (jzh)] J — and a 2 = (X X \ 2 )I, (/A 3 ); J = 0. Here the evaluation of 
two-body t— matrices and permutation operators P, P and P34 as well as coordinate transfor- 
mations due to considering angular momentum quantum numbers instead of angle variables 
lead to more complicated expressions in comparison to our 3D representation. 



V. NUMERICAL TECHNIQUES 

In this section we describe the details of the numerical algorithm for solving the coupled 
F-Y three-dimensional integral equations, Eq. (J5"lj) . The coupled F-Y equations Eq. (j5lj) . 
represent a set of three-dimensional homogenous integral equations, which after discreati- 
zation turns into a huge matrix eigenvalue equation. The dependence on the continuous 
momentum and angle variables (ui,Vi)i = 1,2,3 and x x , x 2 , x X2 , X x , X 2 , X X2 ) is replaced 
in the numerical treatment by a dependence on certain discrete values. Let the numbers 
of these discrete points be denoted by Nj ac , N spi \ and N po i corresponding to momentum 
(ui,Vi)i = 1,2,3), spherical angle (x x , x 2 , X x , X 2 ) and polar angle (x X2 ,X X2 ) variables, the 
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dimension of the eigenvalue problem is: 

N = Nf ac x N% h x N po i x 2 (53) 

The huge matrix eigenvalue equation requires an iterative solution method. We use a 
Lanczos-like scheme, the method of iterated orthonormal vectors (IOV) that is proved to be 



very efficient for nuclear few-body problems 



This technique reduces the dimension of 



the eigenvalue problem to the number of iteration minus one. The eigenvalue equation, Eq. 
f[5"Tj) . schematically can be written as: 

X(E) ip = K(E) i(j (54) 

The kernel of the linear equations K(E) is energy dependent, and \(E) is its eigenvalue with 
if) as the corresponding eigenvector, ip represents the set of F-Y components as ip = (1*). 
For the physical binding energy the eigenvalue X(E) of the matrix kernel K(E) has to be 
one. We start the iteration with two gaussian F-Y components and stop the iteration after 
5-10 times. In order to solve the eigenvalue equation, Eq. fl54l . for the F-Y components, Eq. 
( !5T|) . we use the Gaussian quadrature grid points for the momentum and angle variables. 

The momentum variables have to cover the interval [0, oo]. Because the F-Y components 
drop sufficiently rapidly we limit the intervals to suitable cut-offs. We choose the cut-offs 
as u^ ax = v? ax = v^ ax = v' 3 max , v% ax = u' 2 max and uf ax = v? ax . These cut-off values 
vary depending on the potential we use but they are chosen large enough to achieve cut-off 
independence. 

The iteration of Eq. floTj) requires a three-dimensional and a six-dimensional interpola- 
tion on i(ji(ui U2 U3 x\ X2 x\ 2 ) and a six dimensional interpolation on ^2(^1 ^2 ^3 X\ X2 Xf 2 ) 
for the first F-Y component. Also it requires a five-dimensional interpolation on 
"01 ( u i u 2 u 3 x i x 2 £12) an d a two-dimensional interpolation on ip2(viV2V3XiX 2 Xf 2 ) for the 
second F-Y component. The interpolations should be carry out in the shifted momentum 
and angle arguments. By adding the additional grid points, to all momentum and ±1 to 
all angle grid points, we avoid the extrapolation outside the Gaussian grids. 

We would like to point out that the symmetry properties as shown in Eq. ( 1321) can be 
implemented in the iteration of Eq. ( l5Tl) to cut down the size of the F-Y components and 
thus save time and memory when computing the integrals. 

The functional behavior of K(E) is determined by the two-body t— matrices t s (ui, tv, x; e) 
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and t a (vi, v' 3 , Y 13 >] e*). We solve the Lippmann-Schwinger equation for the fully-off-shell two- 



body t— matrices directly as function of the Jacobi vector variables as described in ref. (42]. 
The so obtained t— matrices are then symmetrized to get t s (ui,u%, x; e) and t s (v%, v' 3 , X; e*), 
where tc, x and Y 13 i in Eq. (lo"Tj) are replaced with new momentum and angle mesh grids Ui,x 
and X. Both angle mesh grids x and X cover interval [—1, +1] and momentum mesh grid 
u\ covers [0,7f ma:c ]. We would like to point out that after having t— matrix t s (ui, u\, x; e) 
(t s (vi, v' 3 , X; e*)) on grids for U\,Ux and x (vi,v' 3 and X) we solve the integral equation 
again to obtain t— matrix at extra points U\ = and x = ±1 {v\ = and X = ±1). 
Thus when iterating Eq. ( |5Tl) we do not have to extrapolate numerically to first momentum 
-Ui(fi) and angle variable x(X) of t s (u\, ui, x; e) (t a (vi,v' 3 , X;e*)). Also we point out that the 
momentum dependencies given in Eq. fl5T|) suggest that we calculate the two-body t— matrix 
x; e) (t s (vx,v' 3 ,X] e*)) for the energies e = E-^-^(e* = E-^-^) dictated 
by the same u<i and u 3 grids [v<i and v 3 grids). So each energy depends on two momentum 
variables. The number of different energies, where both t— matrices are needed, is quadratic 
in the number of momentum grid points. Consequently both t— matrices would be extremely 
huge if we keep the dependence on both momenta. Therefore we introduce two additional 
energy grids e and e* and insert an interpolation step from these grids to e and e*. This 
reduce the memory and computing time necessary for both t— matrices tremendously. In 
order to obtain the second momentum, the angle and the energy for t s (u\,ux,x; e), also the 
angle and the energy for t s (vi,v' 3 , X] e*) required in the iteration of Eq. (IBTl) . we have to 
carry out three- and two-dimensional interpolations respectfully. 

Since the coupled integral equations, Eq. (15T1) . require a very large number of interpola- 
tions, we use the cubic Hermitian splines of ref. 59J for its accuracy and high computational 
speed. It can be useful to mention that in the numerical calculations we use the Lapack li- 
brary 60[], for solving a system of linear equations in calculation of the two-body t— matrices, 
and Arpack library 6l|, for solving the eigenvalue problem. 
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VI. NUMERICAL RESULTS 



A. Three- and Four-Body Binding Energies 



In order to be able to compare our calculations with results obtained by other techniques 
we use the following spin-independent potentials: 



Gauss-type Baker potential 



62] 



V(r) = -51.5 e 



-0.3906 r z 



Gauss- type Volkov potential [63] 

V(r) = 144.86 e~ 1A87r 
separable Yamaguchi potential 

V(p,p r ) 



83.34 e 



[MeV] 



-0.3906 r 2 



[MeV] 



64 1 



— g(p) g(p') ; gip) = 2 \ n 2 



and the spin-averaged Yukawa-type Malfliet-Tjon V potential 65| 



V(r) = 1458.05- 



-3.11r 



-1.55r 



-578.09- 



[MeV] 



(55) 



(56) 



(57) 



(58) 



The parameters used for Yamaguchi potentials are given in table 1. In our calculations with 
above potentials we use m _1 = 41.470 MeV.fm 2 . For four-body(three-body) binding energy 
calculations twenty(thirty two) grid points for angle variables and thirty(forty) grid points 
for Jacobi momentum variables have been used respectively. 

TABLE I: Parameters of the Yamaguchi-type potentials. 



potential 


A[/m" 3 ] 


Plfm- 1 ] 


Y-I 


0.415 


1.45 


Y-II 


0.353 


1.45 


Y-III 


0.182 


1.15 


Y-IV 


0.179 


1.15 



The techniques to which we compare are the VAR 66]- 68] and HEE 69j methods, several 
types of approximating subsystem kernels of the four-body problem by operators of finite 



rank (SKFR) [70|- 



721 ]. the integrodifferential equation approach SIDE and IDEA [73J, the 
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CCE |74j, the ATMS 
I79H and last but not least 2DI 



75|], the GFMC , 



64] 



the DFY 



32|,[77|, the CRCGBV 



the DMC 



80]. 



In table 2 we show the three- and four-body binding energies for Baker potential calculated 
with different methods. Our results for three- and four-body binding energies with values 
—9.76 and —40.0 [MeVj 1 are in good agreement with results of other available calculations. 

TABLE II: Three- and Four-Body binding energies for Baker potential. 
Method 4-body B.E. [MeV] 3-body B.E. [MeV] 



VAR [66] 
VAR [67] 
HHE [69] 
DFY [77] 
DFY [22] 
FY(PW) [35] 



-39.1±0.1 
-40.03 
-40.05 
-40.0 

-39.9989 
-40.03 



-9.76 



FY(3D) 



-40.0 



-9.76 



For Volkov potential our calculations for three- and four-body binding energies yield the 
values —8.43 and —30.2 [MeV] which as shown in table 3 are also in good agreement with 
other calculations. 

TABLE III: Three- and Four-Body binding energies for Volkov potential. 
Method 4-body B.E. [MeV] 3-body B.E. [MeV] 



HH [18] 
SVM [9] 
VAR [67] 
HHE [69J 
DFY [77] 
DFY [32] 
FY(PW) [35 
FY(3D) 



-30.420 
-30.424 
-30.317 
-30.3988 

-30.2 
-30.2467 
-30.27 
-30.2 



-8.43 



-8.43 



The three- and four-body binding energies for separable Yamaguchi type potentials calcu- 
lated with different methods are listed in table 4. Our results for three-body binding energies 
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for Yamaguchi I, II, III and IV are -25.41, -12.45, -9.25, -8.53 [MeV] and for four-body 
binding energies are —89.8, —54.5, —38.2, —36.2 [Mel/], which are in good agreement with 
results of other methods, specially with 2DI. 

TABLE IV: Four-Body binding energies for Yamaguchi type potentials. The numbers in parenthesis 
are three-body binding energies. 

Method Y-I Y-II Y-III Y-IV 



SKFR [70] -84.66 

SKFR [71J -90.10 

SKFR [72] -89.74 
FY(PW) [35] -89.90 (-25.41) 

2DI [64] -89.6 (-25.40) -54.5 (-12.45) -38.3 (-9.24) -36.3 (-8.51) 

FY(3D) -89.8 (-25.41) -54.5 (-12.45) -38.2 (-9.25) -36.2 (-8.53) 



TABLE V: Convergence of the four-body binding energy with increasing number of partial waves 



for Malfliet-Tjon V potential 



4CJ 



h, Ai, A3 


h 


k 


A 2 


Eground[MeV] 














-31.07 


2 


2 








-31.11 


4 


4 








-31.22 


6 


6 








-31.23 


4 


6 


2 





-31.28 


4 


6 


4 





-31.31 


4 


6 


6 





-31.31 


4 


6 


4 


2 


-31.34 


4 


6 


4 


4 


-31.35 


4 


6 


4 


6 


-31.35 


6 


6 


4 


4 


-31.36 


8 


6 


6 


6 


-31.36 



As demonstrated in table 5, the calculation of the four-body binding energy using the 
Malfliet-Tjon V potential in PW scheme converges to value of E = —31.36 [MeV]. Here 
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convergence is reached for li, Ai, A3 = 8 and h,h, A2 = 6 [40|, while the three-body binding 
energy for this potential converges to —7.73 [MeV] 35(. As shown in table 6 our calculations 



for Malfliet-Tjon V yield the va^ 



ue —31.3 [MeV] for four-body binding energy, which is in 



n 



4Qj and SVM |9|] results and with 



good agreement with recent HH [18|, EIHH [29(], F-Y(PW) 
other calculations. Also our result for three-body binding energy with value —7.74 [MeV] 
is in good agreement with the obtained value —7.73 [MeV] of Faddeev calculations in PW 
scheme. 



TABLE VI: Four-Body binding energies for Malfiiet-Tjon V. The numbers in parenthesis are three- 
body binding energies. 



ivietnou 


4-Doay d.fj. [ivievj 


CRCGBV [78] 


-31.357 


ATMS [75] 


-31.36 


GFMC [76] 


-31.3±0.2 


CCE [74] 


-31.24 


VAR [68] 


-31.19±0.05 


IDEA [73] 


-30.98 


DMC [79] 


-31.5 


HH [18] 


-31.347 


SVM [9] 


-31.360 


EIHH [29] 


-31.358 


FY(PW) [35],[4pJ 


-31.36 (-7.73) 


FY(3D) 


-31.3 (-7.74) 



As we can see from these comparisons to other calculations of the four-body binding en- 
ergy based on PW decomposition, our results provide the same accuracy while the numerical 
procedure are actually easier to implement. In the 3D case there is only two coupled three- 
dimensional integral equations to be solved, whereas in the PW case one has two coupled 
sets of finite number of coupled equations with kernels containing relatively complicated 
geometrical expressions. 
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B. Test of Calculations 



In this section we investigate the numerical stability of our algorithm and our 3D repre- 
sentation of Yakubovsky components. We specially investigate the stability of the eigenvalue 
of the Yakubovsky kernel with respect to the number of grid points for Jacobi momenta, 
polar and azimuthal angle variables. We also investigate the quality of our representation 
of the Yakubovsky components and consequently wave function by calculating the expecta- 
tion value of the Hamiltonian operator. For these investigations we use the Malfliet-Tjon V 
potential. 

In table 7 we present the obtained eigenvalue results for binding energy E = —31.3 MeV 
for different grids. We choose the number of grid points for Jacobi momenta as iV Ul = 
N Vl = N V3 = N v > 3 = Nj ac and N U2 = N U3 = N V2 = iV? As demonstrated in this table, 
the calculation of the eigenvalue A convergence to the value one for Nj ac = N 2 ac = 30 
and N sp h = N po i = 20. It should be clear that the solution of coupled Yakubovsky three- 
dimensional integral equations, with six independent variables for the amplitudes, is much 
more time-consuming with respect to the solution of three-dimensional Faddeev integral 
equation [451 ]. with three variables for the amplitude. 



TABLE VII: Stability of the eigenvalue A of Yakubovsky kernel with respect to the number of 
grid points in Jacobi momenta Nj ac , Nj ac , spherical angles N sp h and polar angles N po i, where 
E = -31.3MeV. 



N 1 


N 2 


N sp h = N P oi 


A 


20 


20 


20 


0.987 


26 


20 


20 


0.992 


26 


26 


20 


0.995 


30 


26 


20 


0.998 


30 


30 


20 


1.000 



The solution of coupled Yakubovsky three-dimensional integral equations in momentum 
space allows to estimate numerical errors reliably. With the binding energy E and the 
Yakubovsky components and \1Jj2) available, we are able to calculate the total wave 
function from Eqs. (!33l) - (!35l) by considering the choice of coordinate systems which are 
represented by Eqs. (l39]) - fH0T) . So in order to demonstrate the reliability of our calculations 
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we can evaluate the expectation value of the Hamiltonian operator H and compare this value 
to the previously calculated binding energy of the eigenvalue equation, Eq. fl54l) . Explicitly 
we evaluate the following expression: 

(V\H\V) = (V\H \y) + (y\V\V) 

= ( 12 (V>i|#o|*> + 6<V> 2 |#o|*)) +6(^\V l2 \^f) (59) 



where 



and 



poo poo poo ^2 

(i/j^HqI^) = 8tt 2 / duxul / du 2 u\ \ du 3 ul [— + 
Jo Jo Jo m 



3u 2 2 2uL 

— - H 

Am 3m 



-1 p+1 p2-rr 

x / dx\ I d%2 I dipx ipi(ui u 2 u 3 x\ %2 (fi) ty(ui u 2 u 3 x\ x 2 ip 



1 J -l JO 



poo poo poo 2 2 2 

/+1 r+l r2n 

dX 1 / dX 2 I d(pi v/j 2 (viv 2 v 3 X 1 X 2 (j) 1 ) ^(viv 2 v 3 X 1 X 2 </>i ) 
-1 J-i Jo 



poo /"^l p2iT roo pc 

[^\Vi 2 \^) = 8n 2 / du\ uf / dxi / d(fi / du 2 u\ I dx 2 I 

Jo J-i Jo Jo J-i Jo 

poo 

x ty(ui u 2 u 3 x\ x 2 (fx) / du[uf / dx[ / d(p[ 

Jo J-i Jo 

x V 12 (u 1 ,u' 1 ,y n/ )^(u' 1 u 2 u 3 x' 1 x 2 (p' 1 ) (61) 



du 3 u\ 



(60) 



where yxv = xix[ + a/1 — x\ y/l — cos(</?i — The expectation values of the kinetic 
energy (i^o), the two-body interaction (V) and the hamiltonian operator (H) are given 
in table 8 for Malfliet-Tjon V potential calculated in 3D scheme. In the same table the 
four-body binding energy calculated in 3D scheme is also shown for comparison to the 
expectation values of the Hamiltonian operator. One can see that the energy expectation 
value and eigenvalues E agree with high accuracy. 



VII. SUMMARY AND OUTLOOK 



Instead of solving the coupled F-Y equations in a PW basis, we introduce an alterna- 
tive approach for four-body bound state calculations which implement directly momentum 
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TABLE VIII: The expectation values of the kinetic energy {Ho), the two-body interaction (V) and 
the Hamiltonian operator {H) calculated for Malfliet-Tjon V potential in 3D scheme. 





(H ) [MeV] 


(V) [MeV] 


(H) [MeV] 


E [MeV] 


FY(3D) 


69.7 


-101.0 


-31.3 


-31.3 



vector variables. We formulate the coupled F-Y equations for identical spinless particles 
as function of vector Jacobi momenta, specifically the magnitudes of the momenta and the 
angles between them. We expect that coupled three-dimensional F-Y equations for a bound 
state can be handled in a straightforward and numerically reliable fashion. Our results 
for spin-independent two-body potentials are in good agreement with pervious values for 
VAR, HHE, SKFR and DFY techniques, especially they are matched with PW calculations 
in F-Y scheme. Also working directly with momentum vector variables gives the benefit 
of considering all partial waves, which provides perfect agreement with GFMC, CCE, CR- 
CGBV, ATMS, VAR, IDEA, DMC, HH, SVM, EIHH and F-Y(PW) values for Malfliet-Tjon 
V potential. This is very promising and nourishes our hope that calculations with realistic 
NN potential models, namely considering spin and isospin degrees of freedom, will most 
likely be more easily implemented than the traditional PW-based method. The stability 
of our algorithm and our 3D representation of Yakubovsky components have been achieved 
with the calculation of the eigenvalue of Yakubovsky kernel, where different number of grid 
pints for Jacobi momenta and angle variables have been used. Also we have calculated 
the expectation value of the Hamiltonian operator. This test of calculation has been done 
with Malfliet-Tjon V potential and we have achieved good agreement between the obtained 
eigenvalue energy and expectation value of the Hamiltonian operator. We predict that the 
incorporation of three-body forces will most likely also be less cumbersome in a 3D approach. 
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